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Abstract 

A quasi-geostrophic intermediate complexity model is considered, providing a schematic repre- 
sentation of the baroclinic conversion processes which characterize the physics of the mid-latitudes 
atmospheric circulation. The model is relaxed towards a given latitudinal temperature profile, 
which acts as baroclinic forcing, controlled by a parameter Te determining the forced equator-to- 
pole temperature gradient. As Te increases, a transition takes place from a stationary regime to 
a periodic regime, and eventually to an earth-like chaotic regime where evolution takes place on 
a strange attractor. The dependence of the attractor dimension, metric entropy, and bounding 
box volume in phase space is studied by varying both Te and model resolution. The statistical 
properties of observables having physical relevance, namely the total energy of the system and the 
latitudinally averaged zonal wind, are also examined. It is emphasized that while the attractor's 
properties are quite sensitive to model resolution, the global physical observables depend less 
critically on it. For more detailed physical observables, such as the latitudinal profiles of the zonal 
wind, model resolution again may be critical: the effectiveness of the zonal wind convergence, acting 
as barotropic stabilization of the baroclinic waves, heavily relies on the details of the latitudinal 
structure of the fields. The necessity and complementarity of both the dynamical systems and 
physical approach is underlined. 
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I. INTRODUCTION: ATMOSPHERIC CIRCULATION AS A PROBLEM IN 
PHYSICS AND IN MATHEMATICS 



In the scientific context, Climate is denned by the statistical properties of the Climatic 
System. In its most complete definition, the Climatic system is composed of four intimately 
interconnected sub-systems, Atmosphere, Hydrosphere, Cryosphere, and Biosphere. These 
subsystems interact nonlinearly with each other on various time-space scales jl|, |2 1 . 

The Atmosphere is the most rapid component of the Climatic System. The Atmosphere 
is very rich in microphysical structure and composition and evolves under the action of 
macroscopic driving and modulating agents - solar heating and Earth's rotation and gravi- 
tation, respectively. The Atmospheric Circulation is the basic engine which transforms solar 
heating into the energy of the atmospheric motions determining weather and climate as we 
commonly perceive them. The Atmosphere features both many degrees of freedom, which 
makes it complicated, and nonlinear interactions of several different components coupling 
a vast range of time-space scales, which makes it complex. In many cases, the dynamics 
of such a system is strongly chaotic - in the sense that the autocorrelation function of any 
variable vanishes on finite time scales - and is characterized by a large natural variability on 
different time scales (21 0] 

The understanding of the physical mechanisms operating in the Atmosphere critically 
influences important human activities like weather forecast, territorial planning, etc. This 
is one reason why, more than half a century ago, von Neumann posed the Atmospheric 
Circulation in the core of the ongoing development of numerical modelling jfj. However, 
the General Atmospheric Circulation (GAC) also poses problems of general physical nature 
as a realization - in fact the one we can best observe - of planetary scale thermodynamic 
transformations in a rotating, stratified fluid. 

Historical - see the C .a 8S1 cal olograph and paper by Loren, Q - the problem of 
GAC has been essentially approached in terms of time-mean circu 
which generate and maintain it. Almost one century ago Jeffrey 
to maintain the observed time-mean circulation at middle latitudes, it is necessary to take 
into account the momentum and heat transfer properties of the eddies, i.e. the fluctuating 
component of the atmospheric flows [lo| . 

Among all the physical processes involved in the GAC, the so-called baroclinic conversion 
(baroclinic comes from ancient Greek: constant pressure surfaces not parallel to constant 
density surfaces) plays a central role because it is through this mechanism that rotating, 



ation and the processes 
realized that in order 



stratified fluids convert the available potential energy 



11 



12 



131 ] . stored in the form of 



thermal fluctuations, into the vorticity and kinetic energy of the air flows as we observe 
them. At mid-latitudes of both hemispheres, the baroclinic conversion process can be taken 
as responsible for the destabilization of the fixed point given by the zonally (lon gitu dinally) 
symmetric Atmospheric Circulation characterized by a purely zonal wind (jet) j35|. Baro- 

m 



clinic unstable waves can be actually observed (see for instance 



15, ia). The definition 



of the basic ingredients in the mechanism of baroclinic instability has been one the main 



17 



successes of the dynamical Meteorology of this century 

Within the, virtually innumerable, papers devoted to the subject of GAC, a few happened 
to suggest new methodologies and concepts of general interest for fundamental disciplines, 
such as Physics and Mathematics, as well as more empirical natural and social sciences, such 
as Biology, Medicine and Economics. A leading example is that of Lorenz' attractor jl^j . 
But, apart from such exceptions, the problem of GAC has remained confined within the 
boundaries of Geophysical (mostly Meteorological) literature, with all the ensuing language 
barriers with respect to Physics and Mathematics. Also the relatively recent (last fifteen 
years) public attention on Climate issues has been attracted essentially on phenomenolog- 
ical and/or numerical modelling issues rather than on fundamental mechanisms [^J. One 
consequence of this cultural separation has been that, for example, the knowledge that in 
dynamical systems the stability properties of the time mean state do not even provide a 
zeroth-order approximation of the dynamical properties of the full nonlinear system has 
been, and still is, quite systematically ignored in specialized literature - see |2J| and 



for enlightening examples - despite both theoretical arguments |2Jj] and simple counter 



examples of physical significance 



24L l25| indicated throughout the years. Note that this, 



somewhat methodological, issue bears relevance also in practical problems like the provision 
of the so-called extended range weather forecasts, which extend beyond the deterministic 
predictability horizon of about 10-15 days (see e.g. Lorenz 0, Q). Suppose, in fact, that the 
forecaster was given the next month average atmospheric fields: what practical information 
would he derive from that? Of course, if dynamical information is stored in the average 
fields - for example in the form of dominant regimes of instability derivable from the stabil- 
ity analysis of time-mean flow - we could obtain useful information from the prediction 
of such time mean fields. Unfortunately, as remarked above, this picture is far from being 
true, and the problem of extended range is still open even in terms of formulating clearly 
what we should forecast! 

In order to address some of the above mentioned issues, in this work we consider a 
quasi-geostrophic model of intermediate complexity for the atmospheric circulation. By 

intprmpHifltp wp mpan that thp nnmhpr nf varifl.hlps (AR tn ?>RA\ lips hptwppn thp fpw Hpptppk 



is, tea, and the state-of-the-art Global Circulation 



of freedom of, say, the Lorenz models 

Models j^J, which feature over 10 6 degrees of freedom. The model used here has no seasonal 



cycle and provides an earth-like representation of the turbulent baroclinic jet 2J, |25|. It 



is vertically discretized into two layers, which is the minimum for baroclinic conversion 



to take place 



28l |29| . and latitudinally discretized by a Fourier half-sine pseudo-spectral 



expansion up to order JT. We have used JT = 8, 16, 32, 64, yielding a hierarchy of quasi- 
geostrophic models having increasing resolution. A fundamental property of these models 
is semi-linearity: the eddy field is truncated to one wavenumber in the longitudinal (zonal) 
direction, so that the evolution equation is linear in terms of the time-varying zonal flow. 
This provides a dynamical meaning for the separation between zonal and eddy flow that 
is only geometrical - and originally just geographical - in the traditional approach: in our 
case the zonal flow is an integrator of the nonlinear self-interactions of the wave-field which 
propagates and grows linearly on the zonal flow self. 

In Sec. |n] we present a detailed general derivation of the evolution equations for the 
two-level quasi-geostrophic model starting from the ab-initio equations and explaining the 
approximation involved in the derivation of the 3D quasi-geostrophic equations. This deriva- 
tion allows a clear understanding of the physics involved in the considered hierarchy of 
quasi-geostrophic equations and is alternative to the non-dimensional formulations which 



are common in the meteorological literature 



We further obtain the equations of the 



one-wavenumber model examined in this study, in the form adopted for the numerical inte- 
gration. 

The main results of this work are presented in Sec. IIHI and Sec. IIVI We study the 
sensitivity of the model behavior with respect to the parameter T# determining the forced 
equator-to-pole temperature gradient, which acts as baroclinic forcing. The influence of the 
order of (spectral) discretization in the latitudinal direction is also analyzed. For low values 
of Te there occurs a transition from a stationary to a earth-like chaotic regime. Here chaotic 
means that an attractor is detected having a positive maximal Lyapunov exponent, i.e., a 
strange attractor (see 0] for terminology). In Sec. IIHI we characterize the transition from 
stationary to chaotic dynamics in terms of bifurcation theory and study the dependence 
on Te and on model resolution JT of the dimension of the strange attractor, of the metric 
entropy, and of the volume of its bounding box in the phase space. In Sec. IIVI we analyze 
the statistical properties of two physically meaningful observables, namely the total energy 
of the system and the latitudinally averaged zonal wind. An inspection of the latitudinal 
wind profiles is also presented. In Sec. El we give our conclusive remarks and perspectives 

fnr fntnrp wnrks 



II. THE AB-INITIO FORMULATION OF THE MODEL EQUATIONS OF MO- 
TION 

A. Initial Remarks 

The dynamics and thermodynamics of the dry atmosphere for an observer in the Earth's 
uniformly rotating frame of reference is described by the following equations Q]: 



-^ p + pV-£=0 (1) 
D Vv 

—u + 2Qxu = ^--V$ + F (2) 

P=p(p,T). (4) 



Here p is the density, u is the velocity vector, Q is the Earth's rotation angular velocity, p is 
the pressure, <3> is the geopotential, F is the resultant of the frictional forces per unit mass, 
h is the specific enthalpy, Q is the diabatic heating, D represents the effect of heat diffusion 
processes, and T is the temperature of the fluid. The material derivative D/Dt is defined 
as follows: 

——• = — • + ( -u • V ) • . (5) 



Dt dt 

Equations are commonly referred to as mass continuity, Navier-Stokes, and ther- 

modynamics, respectively. They express the dynamic balances of mass, forces, and specific 
enthalpy of the system, while (j3J) is the equation of state of the fluid under consideration, 
which, in the case of dry air, can be well-represented as a perfect gas. 

The description of the macroscopic behavior of the atmosphere is based on the systematic 
use of dominant balances derived on a phenomenological basis. Suitable approximations to 
equations Q-Q are obtained by assuming that the actual evolution departs only slightly 
from the balances. In fact, different balances have to be applied depending on the time 
and space scales we are focusing on. In this way, it is possible to filter out (exclude) all 
solutions corresponding to physical processes that are heuristically assumed to contribute 
only negligibly to the dynamics of the system, at the time and space scale under examination. 
The magnitudes of various terms the governing equations for a particular type of motion 
are estimated using the so-called scale analysis technique jjyj. The resulting models usually 
five pood annrmnmation tn the ohsprveH fields when sufficiently larpe snatial nr temnoral 



averages are selected 



B. The hydrostatic and quasi-geostrophic approximations 

For the dynamics of the atmosphere at mid-latitudes, on spatial and temporal scales 
comparable with or larger than those of the synoptic weather (about 1000 Km and 1 day, 
respectively), it is phenomenologicall y w ell-established that the hydrostatic balance is obeyed 
with excellent approximation [ll l^H^: 

k-Vp = -pg, (6) 

where: 

k = (7) 

This expresses the balance between the gravitational force and the vertical pressure gradient, 
the vertical direction k being defined by the gradient of <3>. Since the atmosphere is shallow 
with respect to the radius of the Earth, one can use the approximation $ ~ gz, where z is 
the local geometric vertical coordinate. The hydrostatic balance also allows the usage of p 
as vertical coordinate. 

Moreover, in the just mentioned synoptic scales, the atmosphere is close to the geostrophic 
equilibrium, which is realized when the local horizontal pressure gradient exactly balances 
the Coriolis acceleration. In the geostrophically balanced flows, when pressure is taken as 
vertical coordinate, the wind can be expressed as 

u g = (u g , v g , 0) = —k x V$ = A; x V^ 9 , (8) 
Jo 

where /o = 2f2siny? is the the orthogonal projection of the Coriolis parameter on the surface 
of the planet at latitude tp and ip g = $//o is defined as the streamfunction of the flow. The 
geostrophic wind (JHJ) is horizontal and non-divergent. This implies that the geostrophic 
vorticity vector is parallel to the vertical direction and its non-vanishing component can be 
expressed as: 

£ g = k-(V xu g ) = -^-A#$ = A H ^ g , (9) 
Jo 



JflQ 



where Ah is the horizontal Laplacian operator, see e.g. 

Equations © and (JHJ) are only diagnostic, so that no information on the evolution of 
the system can be obtained. From the set of ab-initio dynamic and thermodynamic 



equations of the atmosphere it is possible to obtain a set of simplified prognostic equations 
for the synoptic weather atmospheric fields by assuming that the fluid obeys the hydrostatic 
balance and undergoes small departures from the geostrophic balance. Moreover, we assume 
that the domain is centered at mid-latitudes and it is such that / can be well-approximated 
by the linear expansion / (if) ~ / (<^ ) + 20 cos (<p ) (cp — ip ). 

Local Cartesian coordinates (x, y) and pressure coordinate p are introduced for the hori- 
zontal and vertical directions, respectively, with x denoting the zonal and y the latitudinal 
coordinate. The resulting domain is periodic in x, with wavelength L x , and bounded in y 
and p, yielding 

x G R/2nL x , y G [0, L y ] , p G [0, p ] , (10) 

and / is approximated as / ~ fo + /? (y — L y /2). In the meteorological jargon this is usually 
referred to as the (3- channel. A sketch of the actual geographical area corresponding to the 
/3-channel is presented in Fig. We remark that in this work, in order to avoid problems 
in the definition of the boundary conditions of the system, due to the prescription of the 
interaction with the polar and the equatorial circulations at the northern and southern 



boundary, respectively |24j, we consider a domain extending from the pole to the equator. 
We remark that the quasi-geostrophic approximation is not appropriate for the equatorial 
region, so that we do not expect to capture any realistic feature of the tropical circulation, 
and that the mid-latitude channel is determined by y ranging from 1/4 L y to 3/4 L y , 
corresponding to a latitudinal belt centered about 45°iV with an extension of 45°. 

Proceeding further with simplifying assumptions, the equation of state p = p/RT is 
adopted for (j3J), where R is the gas constant for dry air, so that the following relation holds: 

^ = ~T, (11) 
dp f p 

and the specific enthalpy for the dry air is expressed as h = C p T. We introduce the quasi- 
geostrophic material derivative: 

^=i' + M)"=i" +J <*-">- (12) 

where J is the conventional Jacobian operator defined as J (A, B) = d x Ad y B — d y Ad x B. 
Physically, this means that advection occurs along constant pressure levels and is performed 
by the geostrophic wind. This yields the so-called quasi-geostrophic equations for the stream- 



function ip g : 

(Ai/) g + f + Py)-f -^ = k-VxF + u (A„) 2 ^ (13) 
D g ( dip g \ RpT dQ f Q ( d^ 9 \ R Q 

where u is the velocity in the direction of p, the frictional forces are represented as viscous 
processes with diffusion constant u, the heat diffusion is parameterized by the coefficient k, 
and 6 is the potential temperature: 

R 

Po \ Cp 



Q = T {J) . ( 15 ) 

which is related to the specific entropy s of the air by 

s = C p ln6. (16) 

The quasi-geostrophic approximation is very useful because the resulting evolution equa- 
tions (|13j ) -([14| ) focus on the process of slanted convection which is responsible both for the 
baroclinic conversion of potential energy into eddy energy and for the generation of vorticity. 
These are the essential ingredients underlying the generation of atmospheric disturbances 
at mid-latitudes [l|, 0, Q, Q 

The non-geostrophic velocity component u does not have an evolution equation and can 
be diagnosed from the thermodynamics equation (|14j). The following boundary conditions 
apply for the ageostrophic velocity uj: 

u{x, y,p = 0) = (17) 
w(x,y,p = po) = -E £g(x,y,p = po) (18) 



where the condition at p = po is due to the E 
atmosphere with the planetary boundary layer 
approximation: 



;man description of the coupling of the free 
. We adopt the phenomenologically-based 



^HZ - -H(P)\ (19) 



RpTdB 

where H (p) is the vertical scale related to the stratification of the atmosphere which de- 
pends only on p. By substituting the vertical velocity u obtained in equation (|14|) into 



equation (fT^j) . we obtain that the quasi-geostrophic potential vorticity q g , defined as: 

„ = + + (20) 

satisfies the following canonical equation 



32. 



3: 



Dt 



K- 



d_ 

dp 



p 



dp 







(21) 



The quantity q g (as well as all of its powers) is conserved along motion if no diabatic forcing 
is applied (Q = 0) and if the diffusion and viscous effects are discarded, e.g. by setting in 
our case v = k = 0. 

We remind that, formally, the quasi-geostrophic equations ()13|) -([14j ) can be derived from 
the ab-initio equations (HJ-(HJ) by retaining the zeroth and first order term in the expansion 
performed on the Rossby number: 



Ro 



U 



(22) 



with (3L <ti /q, where U and L are typical values of the horizontal velocity and horizontal 
space scale [28]. 



C. The two-level model 

A simplified version of system (|T3j) and ((H)) is produced by discretizing the vertical di- 
rection into a finite number of pressure levels. This vertical discretization approach has been 
first introduced by Phillips 29( and retains the baroclinic conversion process, which is the 
basic physical feature of the quasi-geostrophic approximation. 

We refer to Fig. El for a sketch of the vertical geometry of the two-layer system. The 
streamf unction if} g is thus defined at pressure levels p = pi = Po/4 and p = pz = 3/4p , while 
lj is defined at the pressure levels p = (top boundary), p = p 2 = po/2, and p = po (surface 
boundary). The pressure level pertaining to the vertical derivative of the streamfunction 
dipg/dp as well as the stratification height H is p = p 2 - We note that 5p = p% — p\ = p 2 = 



Po/2- This system is described by the following equations of motion: 



^ (A H fa + fo + 0y) - fo^jf^ = 0, (23) 
^ (A^ 3 + fo + fa) - fo^J^ 1 = 0, (24) 



D 2 ffa-fa\ 2 f (fa-fa\ R Q2 , _x 

where we have neglected the viscous dissipation by setting v = 0, dropped the g subscript 
for simplicity, and have adopted the notation 

fa = J = 1,3, Q 2 = Q(p 2 ), (26) 

u 3 = u {pj) , j = 0, 2, 4, tf 2 = #(p 2 ), (27) 

S # = Jr +J( ^' #) J = 1 ' 3 - (28) 
The boundary conditions ([17| ) -([18j ) on the vertical velocity are implemented as 

^0 = 0, (29) 
u 4 = -E A H fa, (30) 

where the streamfunction at the top of the boundary layer has been approximated by the 
streamfunction ^3 



28]. The streamfunction at the intermediate level p 2 is computed as 
average between the streamfunctions of the levels p\ and p^, so that the material derivative 
at the level p 2 can be expressed as: 



Dt 2\Dt Dt) dt 2 v</y1 ^^ 3 ' > v ' 



Along the lines of the derivation of (|21| ) -()20 p . by substituting u 2 , u , and uj^ (as defined in 
(|23j). (|2l?|). and (J3UJ), respectively) into and (jUJ), we obtain the evolution equations for 
the quasi-geostrophic potential vorticity at the two levels: 



Dl K a ! 1 1 \ R Q2 t ~ , 

Di qi = -Hl Ah ^ - ^ ~ J^-C^ (32) 



Here the q^s are defined as: 



q t = A^i + fo + Py + (1 - 2^,0 -2 (Vi - Va) , i = 1, 3, (34) 

where 5^1 is the Kronecker's delta, which is equal to 1 when the two indexes are mutually 
equal and otherwise. 

It is possible to derive the following expression for the horizontal energy density of the 
system: 



e (x, y) = — 
9 



(35) 



[ 2 

Here the factor Sp/g is the mass per unit surface in each level, the last term and the first 
two terms inside the brackets represent the potential and kinetic energy, respectively, thus 
featuring a clear similarity with the functional form of the energy of a harmonic oscillator. 
We emphasize that in (f3*H|) the potential energy term is half of what reported in [3], which 
contains a trivial algebraic mistake in the derivation of the energy density, as discussed with 
the author of the book. 

We choose the following simple functional form for the diabatic heating: 

Q 2 = VN C p{r -T^ VN C p ff(^-^.), (36) 

where r* has been introduced for later convenience and, consistently with equation (jllj) . T2 
is evaluated at the pressure level 2 and is defined by 

^3-^i = _J_t 2 . (37) 

<5p f P2 2 ' 

The functional form of equation (|36J) implies that the system is relaxed towards a prescribed 
temperature profile T* with a characteristic time scale of 1/Vjv- T* and r* are respectively 
defined as follows: 

rp* t e f ny\ * RT E { ny\ {QQ , 

T =- C08 U)- T =AT C08 Uj- (38) 

so that Te is the forced temperature difference between the low and the high latitude border 
of the domain. Since we assume no time dependence for the forcing parameter Te, we 
discard the seasonal effects. Considering that the thermal wind relation: 

^=t x a (39) 

op op 



can be discretized as follows for the two level system: 

= k x V^f^, (40) 
op Op 

we have that the diabatic forcing Q2 in (}3T)j) causes a relaxation of the vertical gradient of 
the zonal wind u\ — u 3 towards the following prescribed profile 2m*: 

R 7T T E . / 7iy\ 
2m = ———sin — , (41) 



fo L y 2 \ L 

where the constant 2 has been introduced for later convenience. We introduce the baroclinic 
and barotropic components (r, 0) as 

r=\{?k- *h) , (42) 
0=^(^1 + ^3). (43) 

From equations (|32 ft - (j33j) one obtains the equations of motion for (r, </>): 
d 2 8/ 2 \ 

dt AHT ~lPdt r + J \' Ah<I) + f3y + IP v + J ^ A//r) = 

/ , s 2k a 2z/jv . ,. . A A . 

-j^A H (0 - r) - ^A ff r + (r - r*) , (44) 



^ A , T / I A / „ N X / A \ 2l/ E 



^-A H + J (<p, A H <P + /3y) + J (r, A H r) = ~^&H (<P ~ r) . (45) 



where ^ = foEoH%/ (26) and the meaning of r* is made clear. Notice that this system only 
features quadratic nonlinearities. The two- level quasi-geostrophic system (jS|)-((13J) can be 
brought to the non-dimensional form, which is more usual in the meteorological literature 
and is easily implementable in computer codes. This is achieved by introducing length 
and velocity scales I and u and performing a non-dimensionalization of both the system 
variables (x,y,t,<fr,T,T) (as described in Table HJ) and of the system constants (Table [HJ). 
When assessing, as in our case, atmospheric phenomena from synoptic to planetary scales, 
suitable choices for the length and velocity scales are I = 10m 6 and u = 10ms -1 . With the 
choices of the constants described in table |H1 our system is equivalent to that of Malguzzi 



and Speranza 



24j . where the following correspondences hold: 



1 „ 2is E v E 2k 
Jp^ F ' ~H2**~2' Hl^ USl VN VH ' ^ ' 



D. The single zonal wave two-level model 



In this section we derive the evolution equations used in the present study. We Fourier- 
expand the and r fields in the zonal direction x as follows: 

oo 

<fr (x, y,t) = ^ A n (y, t) exp (i2nnx/ L x ) + c.c. (47) 

n=0 

oo 

r(x,y,t) = ^ B n (y, t) exp (i2mrx/L x ) + c.c, (48) 

n=0 

where c.c. stands for complex conjugate. By definition we have 

U(y,t) = -^!M., m (y,t) = - a -*Ml, (49) 

so that U represents the zonal average of the mean of the zonal wind at the two pressure levels 
1 and 3 (see previous Section), while m represents the zonal average of the halved difference 
between the the zonal wind at the two pressure levels 1 and 3. In this work we focus on 
the interaction between the average zonal wind and waves, thus neglecting the wave-wave 
nonlinear interactions. We therefore only retain the zonally symmetric component (i.e., that 
of order n — 0) and one of the non- zonal components (i.e., for a fixed n > 1) in the Fourier 
expansions (|47|) -(|48| ) and in the equations of motion. Since quadratic nonlinearities like those 
described in equations ()44j) - (|43|) generate terms with Fourier components corresponding 
to the sum and difference of the Fourier components of the two factors, no wave-wave 



interactions can take place. Note that i 
interaction would have been possible 



cubic nonlinearities were present, direct wave- wave 
. In the present case, the wave can self-interact only 
indirectly through the changes in the values of the zonally symmetric fields U and m. This 
amounts to building up a semi-linear equation for the wave on top of a nonlinear dynamics 
for the zonally symmetric parts of the fields. 

As the only retained non-zonal component we select that of order n = 6, since we intend 
to represent the baroclinic conversion processes, that in the real atmosphere take place on 
scales of L x /6 or smaller Q|. With this choice, setting \ = 6 x 2tt/L x in order to simplify 
the notation, equations (JH|)-(jlBI) reduce to 

ry 

<p(x,y,t) = — I U (z, t) dz + Aexp (ivx) + c.c, (50) 

Jtt/2 

ry 

r(x,y,t) = — I m (z, t) dz + 5 exp (ivx) + c.c, (51) 

Jn/2 



where the choice of the lower integration limit will be explained later. By substitut- 
ing (p^Ujl - lfoTj) into equations (jHJ)-(jlHl) and projecting onto the Fourier modes of order n = 
and n = 6, we obtain the equations: 



(52) 



Ayy ~ + (i X U + Ayy ~ ^ U + ^Uyy + ^ ^ ~ i)(P\ A 

( 2ve\ / n 2ve o\ n 
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where (I52j) - (J53|) and (|54j) - (|55j) refer respectively to the non-zonal and zonal components, the 
dot indicates time differentiation, and X* denotes the complex conjugate of X. This is a set 
of 6 equations for the real fields A 1 , A 2 , B 1 , B 2 , U, m, where A 1 and A 2 are the real and 
imaginary parts of A and similarly for B. Rigid walls are taken as boundaries at y — 0, L y , 
so that all fields have vanishing boundary conditions. We emphasize that, by construction, 
no wave- wave interactions occur in ()52jl - ()55jl . Moreover, only quadratic nonlinear terms are 
present, due to the fact that the same holds for (|44j) -(|45| ) . A Fourier half-sine expansion of 
the fields is carried out, with time- varying coefficients: 
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E m ^' sin ( 1~) > ( 59 ) 



truncating at order JT. Therefore, the lower integration limit in (|5() j) - (j51|) is such that 



the two fields r and <fi as reconstructed from (|56 p -([59 p have automatically zero mean when 
latitudinally integrated. Such a choice allows for the fact that the energy density (|35|). and 
consequently the total energy of the system, does not depend on the latitudinally averaged 
value of r, which has no physical relevance. We denote by IX, (•) the projection operator 
onto the basis function sin(^). Because of computational speed, we chose a collocation 
(also known as pseudospectral) projection, see Appendix [X] for a details. By linearity of 
IT, (■), its action on linear terms in (|52j ) -(|55 p is obvious. For example, terms like A y are 
represented as 

A«/ = — / w 2 A) sin I — — ] , where Wj = — . (60) 

W Pi 3 3 \L„J Ly 

So by plugging expansion (|55 j) -(|55 )) into the equations (|52 |l -(|55 jl . and by applying ILj (■), we 
eventually obtain a set of 6 x JT ordinary differential equations in the coefficients Aj, A?, 
Bj, B 2 , Uj, rrij, with j = 1, . . . , JT: 
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System (|6T ]) -(f6T H) constitutes the base model of our study. For the truncation order JT we 
have used the values: JT = 8, 16, 32, 64. 



III. DYNAMICAL AND STATISTICAL CHARACTERIZATION OF THE 
MODEL'S ATTRACTOR 

A. Hadley Equilibrium 

The system of equations (|44 |) -(|45 |) has the following stationary solution for zonally sym- 
metric flows: 



<P(y)=r(y), (67) 
2k d 2 r (y) i 2u N 

l 2 



Considering the functional form (|38|) for r* (y), the following expression for r (y) holds: 



RT E c ™{- y T *(y) 

r(y) = <P(y) = T - A ~r\2= TT*- (69) 

Jo 4 i + (jl) l + JL) 

U N \Ly J U N \Ly J 

When expressing this solution in terms of the average zonal wind U (y) and of half of the 
wind shear m (y), both defined in (|49|). we have: 

m(y) , U(y) K^J^ML = ^M- i , (TO, 

UN \Ly J U N \Ly J 

where we have used the definition ([41 p for m*(y). Moreover, since the temperature T 2 
is proportional to r (compare (|3*7j) and (j42| ). the following temperature profile T 2 (y) is 
realized: 

T 2 (y) = —-a = —-2- (71) 
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This solution describes a zonally symmetric circulation characterized by the instantaneous 
balance between the horizontal temperature gradient and the vertical wind shear, which 
corresponds on the Earth system to the idealized pattern of the Hadley equilibrium 28, 
351 ] . In particular, since m (y) = U(y), we have that u\ (y) = 2U (y) and u 3 (y) = 0, 
i.e. all the dynamics takes place in the upper pressure level. Since the the lower pressure 
level experiences no motion, the Ekman sucking process is switched off and, consistently, 
the solution does not depend on the corresponding coupling constant v E - 

There is a value of the equator-to-pole temperature gradient Tj? such that if Te < T^ 



the Hadley equilibrium (|67 p -([70 p is stable and has an infinite basin of attraction, whereas if 
Te > Tg it is unstable. 

In the stable regime with Te < T^ , after the decay of transients, the fields 0, r, m, U, 
and T are time-independent and feature zonal symmetry - they only depend on the variable 
y. Moreover, they are proportional by the same near-to-unity factor to the corresponding 
relaxation profiles, compare (|67j ) -([71| ) . In particular, this implies that all the equilibrium 
fields are proportional to the parameter Te- 

In our model, since the forcing m* (y) only projects onto the first latitudinal Fourier mode 
(see (JHJ)), the Hadley equilibrium is fully described as follows: 

A){t)=B i j (t)=O t z=l,2, j = l,...,JT, (72) 
m j (t) = U j (t)=0, j = 2,...,JT, (73) 

mi (y) = Ut (y) = m * (y) x2 - (74) 
1 + — ( f- 

When increasing the values of the control parameter Te beyond , the equilibrium 
described by (J57j) -(|55 |) becomes unstable. The physical reason for this, as first pointed out 
by Charney and Eady on vertically continuous models Q] and by Phillips on the two 
two- level model 3], is that for high values of the meridional temperature gradient the 
Hadley equilibrium is unstable with respect to the process of baroclinic conversion, which 
allows the transfer of available potential energy of the zonal flow stored into the meridional 
temperature gradient into energy of the eddies, essentially transferring energy from the latter 
term to the first two terms of the energy density expression ([35)1 . The two-level model, as 
first pointed out by Phillips 2^1, is the minimal model allowing for the representation of 
this process. 

At mathematical level, in our model we have that for Te = T^ a complex conjugate pair 
of eigenvalues of the linearization of (|61| ) -(|66 j) cross the imaginary axis so that their real part 
turns positive, which suggests the occurrence of a Hopf bifurcation j^fj. 

The observed value of Tjf changes with the considered truncation order JT. Results are 
reported in Table ITTT1 for the choice of constants reported in Table ITT1 and for JT = 8, 16, 



32, 64. We have that T% increases with the value of JT. The reason for this is that the 
finer is the resolution, the more efficient are the stabilizing mechanisms which counteract the 



baroclinic instability. Such mechanisms are the barotropic stabilization of the 



the horizontal shear through the convergence of zonal momentum 
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et, increasing 



39,Ufl, and the 



viscous dissipation, which both act preferentially on the small scales since they involve the 



spatial derivatives of the fields <fi and r. This is a clarifying example that in principle 
it is necessary to include suitable renormalizations in the parameters of a model when 
changing the resolution JT, in order to keep correspondence with the resulting dynamics j^. 
Nevertheless, in our case the values of obtained for the adopted resolutions are rather 
similar. 

Moreover, in our model the number of linearly unstable modes of the Hadley equilib- 
rium ([67 )1 - ([70j ) increases with the value of T#. As shown in Fig. El at each jump in the 
graphs an additional pair of complex conjugate eigenvalues crosses the imaginary axis. The 
increase of the number of linearly unstable modes of the Hadley equilibrium can be framed 
at physical level in the fact that for larger values of Te a larger pool of available potential en- 
ergy is available for conversion and faster latitudinally varying modes can become unstable, 
similarly to case of the Phillips model j^J . We remind that the system under investigation 
obeys a Squires condition [41 1 . so that the fastest growing among the unstable modes is the 
latitudinally gravest one. The signature of the relevance of the stabilizing mechanisms and 
of the geometrical properties of the linearly unstable modes developing for higher values of 
JT can be confirmed by observing respectively that while for low values of Te (Te ^$ 12) the 
number of linearly unstable modes decreases with JT, the converse is true for high values 
oiT E (T E >25). 



B. Transition to Chaos 



In this section we analyze the route to the formation of a strange attractor of 
model (jo f T|) -(jol) j) as the parameter Te is increased. Throughout the section, JT is fixed 
at 32, but the results are similar for the other considered values of JT. 

A stable periodic orbit (Fig. 0] (A)) branches off from the Hadley equilibrium (jri7j) -(|7() jl 
as Te increases above Tjf ~ 8.275. We recall that the Hadley equilibrium loses stability at 
Tj? as a pair of complex conjugate eigenvalues of the linearization of (JfTTf - lJBl)!) crosses the 
imaginary axis, (see previous section). This strongly suggests the occurrence of a super- 
critical Hopf bifurcation, which might be checked by center manifold reduction and normal 
form analysis (see e.g. [36]), but it is beyond the scope of the present work. The attracting 
periodic orbit persists for Te in a narrow interval, up to approximately Te = 8.485, where 
it disappears through a saddle-node bifurcation takin g pl ace on an attracting invariant two- 
torus, see Fig. Intermittency of saddle-node type |42J on the two-torus is illustrated in 
Fig. |3](B): after an initial transient the orbit is attracted to a quasi-periodic evolution char- 
acterized by long time-spans, resembling the periodic evolution of Fig. El (A), alternated by 



relatively short bursts in which the orbit explores the rest of the two-torus. In other words, 
for Te right after the saddle-node bifurcation, the orbit on the two-torus slows down in the 
phase space region where the saddle-node has taken place. This yields a higher density of 
points in that region, see Fig. El 

For slightly larger values of Te, a strange attractor develops by so-called quasi-periodic 
breakdown of a doubled torus. This is one of the most typical routes for onset of chaos 



43. 



(weak turbulence) in fluid dynamics experiments and low-dimensional models, compare 

and references therein. We describe this route by means of a Poincare section of 
the attractor of (pT|) - (ffifi|l . obtained by intersecting an orbit with a hyperplane U± — cq for 
a suitable constant Cq. In this Poincare section, the two-torus yields a circle (Te = 8.516) 
which is invariant and attracting under the Poincare (return) map, see Fig. El (A). At first, 
at Te = 8.52 the two-torus loses stability through a quasi-periodic period doubling (see j^, 
Sec. 4.3] and references therein for the theory of quasi-periodic bifurcations). Thereby a 
period two circle attractor is created, meaning a pair of disjoint circles mapped onto each 
other by the Poincare map (Fig. El (B)). By further increasing Te up to Te = 8.521, a 
second doubling occurs, where a period four circle attractor is born (Fig. El (C)). Then for 
Te = 8.522 approximately a transition to chaotic motion occurs: the period four circle turns 
into a strange attractor having a narrow band-like structure (Fig. El (D)), which is likely to 
be a quasi-periodic Henon-like strange attractor, see [44], We remark that: 

• For smaller values of Te, a sort of doubling bubble occurs, i.e. two consecutive dou- 
blings (at approximately Te = 8.504 and Te = 8.509) resulting in an attractor like in 
Fig. El (C), followed by two undoublings (at T E = 8.513 and T E = 8.516) where the 
circle attractor in Fig. El (A) reappears. This sort of direct-inverse finite sequence is 
not uncommon in dynamical systems, see e.g. j^ . 

• Both the breakdown of a quasi-periodic circle attractor and the resulting quasi-periodic 
Henon-like strange attractor are dynamical phenomena occurring rather frequently but 
which are not completely understood from the theoretical viewpoint (see 44j,|48j]). 



As T E further increases, the band widens and blurs (Fig. El (E), for T E = 8.58) until no 
significant structure can be visually detected (Fig. El (F)) for Te = 10. 

The statistical properties of the attractor of (|6i p - (j66|) also display a typical evolution. 
For values of Te nearby the two-torus breakdown, quasi-periodic intermittency is observed, 
i. e. the autocorrelations of an observable (a function of state space variables) typically decay 



very slowly. We consider the total energy E(t) of (|6i p -([66 p . defined as: 



E(t) — / / e(x,y,t)dxdy = 6 / / e(x,y,t)dxdy, (75) 
Jo Jo Jo Jo 

where e(x,y,t) is the energy density in (jSSJ) (details on the algorithm used for the compu- 
tation of time series of E(t) are given in Appendix EJ). In Fig. (A), we display the lagged 
autocorrelation ACF[E(t), Lag] of the time series of the total energy for Te = 8.521, 8.522 
and Te = 8.58. These cases are representative of the qualitatively distinct observed behav- 
iors. For Te = 8.522 the autocorrelation is very similar to what obtained for Te = 8.521, in 
spite of the fact that the former value corresponds to chaotic behavior whereas the latter to 
regular (quasi-periodic) dynamics. This occurs because for Te = 8.522 the chaoticity is very 
weak and the quasi-periodic intermittency rather strong. A much faster decay of the auto- 
correlation, albeit still with the signature of intermittency, is observed for Te = 8.58. The 
quasi-periodic intermittency for Te near T E rit is also illustrated by the geometrical structure 
of the attractor, which still bears resemblance with that of the formerly existing torus. See 
the Poincare sections in Fig. El (D) and (E). 

For larger values of Te (Fig. Ej), we have that at Te = 9 the autocorrelation decays quite 
similarly to the case Te = 8.58 (the Poincare section, not shown, is also quite similar), 
but already at Te = 10 the quasi-periodic intermittency is no longer present (compare 
Fig. El (F)). Correspondingly, the autocorrelation dec ays rather quickly for Te = 10 and, a 
fortiori, for Te = 18. Again compare with j3, 44 , 4^1 [ 
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Completely analogous routes to chaos occur for model (|6T|) - (f6T)j) with .IT = 16 and 64. 
However, the locations on the T^-axis of the various bifurcations are slightly shifted with 
respect to the case JT = 32, compare Fig. Eland Table UTTl Moreover, for JT = 8, a different 
route takes place, involving a quasi-periodic Hopf bifurcation of the two-torus (instead of a 
quasi-periodic period doubling), whereby an invariant three-torus is created. In the Poincare 
section (not shown), this corresponds to an attracting two-torus. 

The invariant objects involved in the transition to low-dynamical chaos described in this 
section correspond to well-known fluid flow patterns. In particular, the two-torus attractor 
in phase space yields an amplitude vacillation in the flow, whereas the three-torus detected 



for JT = 8 yields a modulated amplitude vacillation, see |46|, |50| and references therein. 



However, a characterization of the strange attractors occurring for large Te and of their 



relation to turbulence is still lacking. Typically, low- dimensional nonhyperbolic strange 
attractors, such as the Lorenz [19] and Henon-Pomeau attractors are the topological 

closure of a set of unstable periodic orbits. Moreover, the Henon-Pomeau attractor coincides 



with the closure of^the unstable manifold W u (p) of a fixed point p of saddle type. See 
e.g. 



[ y, y, 

references therein. To the best knowledge of the authors, no 
similar properties has yet been proved (or even formulated) with sufficient generality for 
nonhyperbolic strange attractors of larger dimension. 

We suspect that plenty of unstable periodic orbits and invariant tori coexist with the 
attractor of model (|61 )) -([66 )) with JT = 32, for sufficiently large Te- Indeed, from Fig. |3]we 
deduce that the Hadley equilibrium undergoes several other bifurcations after the first one. 
Since the number of unstable eigenvalues of the Hadley equilibrium increases at each Hopf 
bifurcation, the periodic orbits that branch off have unstable manifolds of increasingly high 
dimension. Moreover, these unstable periodic orbits in turn undergo Hopf bifurcations (also 

n n 

called torus or Neimark-S acker |36j) where unstable two-tori branch off, compare |2J, Sec. 
5]. It seems, therefore, that the phase space quickly gets crowded with high-dimensional 
unstable invariant manifolds. The question remains open whether such complex dynami- 
cal characterizations of the system play a role in the geometrical structure of the strange 
attractor and are potentially useful for computing the statistical properties, let it go for 
the time average fields considered in the classical atmospheric circulations theories or the 
Hadley equilibrium of most theories of atmospheric instability. 



C. Lyapunov Exponents and Dimension of the Strange Attractor 

To characterize the dynamical properties of the strange attractors of (|61 |) -()66 |) we resort 
to the study of the Lyapunov exponents 0, 0] . See Appendix |B] for a description of the 
algorithm used to compute them. In what follows, the Lyapunov exponents are denoted by 
Ai, A 2 , . . . , \ N , with Ax > A 2 > ■ ■ ■ > X N , N = 6 x JT. 

In the left panel of Fig. [HI we represent the evolution of some of the 192 Lyapunov expo- 
nents of the attractor of ()oT)) - ()fifij) with JT = 32 as Te is increased. The maximal exponent 
Ai becomes positive as Te crosses the torus breakdown value T E rit , and then increases mono- 
tonically with T#. 

The spectrum of the Lyapunov exponents is plotted in the right panel of Fig. |S]for three 
different values of Te, again with JT = 32. The distribution of the exponents approaches a 
smooth shape for large Te and a similar shape is observed for JT = 64 (not shown). This 
suggests the existence of a well-defined infinite baroclinicity model obtained from (JUJ-tJUiJ) 
as a (possibly, singular perturbation) limit for Te —> oo. 



1. Dimension of the Strange Attractor 



The Lyapunov exponents are used to compute the Lyapunov dimension (also called 
Kaplan- Yorke dimension, see 13,157 1) and metric entropy (also known as Kolmogorov-Sinai 
entropy j^])- 

The Lyapunov dimension is defined by 

D L = k + ^f^, (76) 
|^fc+i| 

where k is the unique index such that $2j=i — an d 2^=i < 0- Under general 
assumptions on the dynamical system under examination, Dl is an upper bound for the 
Hausdorff dimension of an attractor. 

We have also computed (not shown) other numerical estimates for the dimension of an at- 
tractor: the correlation and information dimensions However, these estimates become 
completely meaningless when the Lyapunov dimension increases beyond, say, 20. In partic- 
ular, the correlation and information algorithms drastically underestimate the dimension. 
This is a well-known problem: for large dimensions, prohibitively long time series have to be 
used 

0- Quelle Q suggests the following rule of thumb: you need a time series of length 
10 d / 2 to estimate an attractor of dimension d. Therefore, computational time and memory 
constraints in fact limit the applicability of correlation-like algorithms to low-dimensional 
attractors. 

The number of positive Lyapunov exponents (unstable dimension JscJ) increases with Te, 
which implies that the Lyapunov dimension also does so. This is confirmed by a plot of the 
Lyapunov dimension as a function of Te for four values of the discretization order JT = 8, 
16, 32, and 64 (see Fig. EJ). For all the considered values of JT, it is possible to distinguish 
three characteristic regimes in the behavior of the function Dl(Te)' 

• For small values of (Te — T E r%t ), we have that Dl oc (Te — T|T**) 7 , with 7 ranging from 
~ 0.5 (JT = 8) to ~ 0.7 (JT = 64). The range of Te where this behavior can be 
detected increases with JT. 

• For larger values of Te a linear scaling regime of Dl ~ (3Te + const, is found in all 
cases. The linear coefficient is for all JT remarkably close to (3 ~ 1.2. The domain of 
validity of the linear approximation is apparently homothetic, as can be seen from the 
simple geometric construction in figure Fig. El 

TV, 1 q r o or TT.A 



saturation as the Lyapunov dimension begins to increase sublinearly with Tg. Note 
that while for JT = 8 the model is in this regime in most of the explored T^-domain 
(Te > 20), for JT = 64 the threshold is reached only for Te > 108. In this latter regime 
of parametric dependence the system is not able to provide an adequate representation 
of the details of the dynamics of the system. Further discussions on this point will be 
given in Sec. IHI Dl and Sec. IIVI 



2. Entropy production 

The metric entropy h(p) of an ergodic invariant measure p expresses the mean rate of 



information creation, see 



for definition and other properties. If a dynamical system 



possesses a SRB ( Sinai- Ruelle-Bowen) invariant measure p, then Pesin's identity holds: 

Hp) = J2 V ( 77 ) 

A,->0 

Existence of an SRB measure for is rather difficult to show for a given nonhyperbolic at- 
tractor 30]. It has been only proven for low-dimensional cases such as the Henon 0| or 
Lorenz |55( strange attractors. More generally one has the inequality h(p) > J2 Xj >o^r We 
then simply assume the existence of a unique SRB measure and refer to the sum of the 
positive Lyapunov exponents as metric entropy. 

The maximal Lyapunov exponent, the predictability time t p = X^ 1 , and the metric en- 
tropy as functions of Te are compared for JT = 8, 16, 32, and 64 in Fig. It turns 
out that, for fixed JT, Ai increases sublinearly with Te, whereas for Te fixed, Ai decreases 
for increasing values of JT. Consequently, for fixed JT the predictability time decreases 
monotonically with Te- We note that, for all values of JT, if Te > 14 we have that t p < 10, 
which corresponds in physical units to a predictability time t p < 12 days. Moreover, in the 
range T E > 12, t p is proportional to {T E — T|P*) 7 , with 7 ranging between [—0.85, —0.8] de- 
pending on the considered value of JT. The metric entropy has a marked linear dependence 
h ~ (3{T E - T E nt ), with (3 ranging from ~ 0.15 (JT = 8) to ~ 0.5 (JT = 64). Moreover, for 
a given value of Te, the metric entropy increases with JT. From the dynamical viewpoint, 
this means on one hand that the maximal sensitivity of the system to variations in the initial 
condition along a single direction is largest for JT = 8. On the other hand, there are many 
more active degrees of freedom for JT = 64 and they collectively produce a faster forgetting 
of the initial condition as time goes on. 



3. Parametric smoothness of the attractor properties with respect to Te 



The dependence of the Lyapunov exponents and, consequently, of the predictability time, 
of the Lyapunov dimension and metric entropy, with respect to T# is remarkably smooth, 
especially if one keeps in mind the paradigms of low- dimensional nonhyperbolic strange 



attractors. For example, for the logistic mapping (see e.g. |3Cj) the maximal Lyapunov 
exponent Ai is a discontinuous function of the parameter at every point where Ai > 0. This 
is due to the fact that so-called windows of periodicity, that is, open parameter intervals 
where the logistic mapping has a periodic attractor, are dense in the parameter axis. In 
the complement set of the windows of periodicity, parameter values for which a strange 
attractor occurs form a nowhere dense set of positive Lebesgue measure. In fact, similar 
features seem to hold for many low- dimensional ma ppin gs having strange attractors, such 



as the Henon-like families 



52 



n 



and references therein. 



54 |6Cj , also compare [30|, |4j 
No windows of periodicity were detected in the fully chaotic range (say, Te > 16) for 
model (j6Tj) - (j66j) . independently of the truncation order JT = 8, 16, 32, 64. We have also tried 
slightly different spectral discretization schemes and integration methods (such as leapfrog 
or Runge-Kutta 4), but this qualitative feature of smoothness and absence of windows of 
periodicity persisted in all cases. 

There are two possible explanations for this: either the windows of periodicity are very 
narrow or there are no windows of periodicity. A possible theoretical support for the latter 
case might be provided by the concept of robust strange attractors. We refer the interested 



reader to 



6l| 



for a class of low-dimensional 



55| for a discussion and more references. Also see 
maps where strange attractors occur on open parameter sets. 

From the above it follows that, from the dynamical point of view, the model (p)Tf - (j6T>l) 
behaves in sensibly different ways if the truncation order JT is changed. For example, in 
the earth- like regime Te = 18, the Lyapunov dimension nearly doubles when passing from 
JT = 32 to JT = 64. However, despite the quantitative differences, many qualitative 
features remain the same for JT > 16: 

• the route for the creation of the strange attractor involves a Hopf bifurcation of the 
Hadley equilibrium, followed by quasi-periodic breakdown of the invariant torus; 

• a linear scaling regime exists for the Lyapunov dimension as a function of Te] 

• the maximal Lyapunov exponent and the metric entropy increase monotonically with 



the distribution of the Lyapunov exponents tends to a well-defined shape for Te large 
(Fig. H right); 



• the dependence of Lyapunov exponents, dimension and metric entropy with respect 
to Te is remarkably smooth. 



D. Bounding Box of the Attractor 

In this section we study the volume of the bounding box Vbb for the attractors 
of model (j61|) - (|UU|) previously described. The bounding box of a set of points in an 
N— dimensional space is defined as the smallest hyperparallelepiped containing the con- 
sidered set 0,0]. For clarity, in the ^-dimensional phase space, where N = 6 x JT, the 
volume Vbb is computed as: 



N=e xJT 



V BB = 



k=l 



max (z k (t)) - min (z k (t)) 



(78) 



Here the z k denote the 6 x JT variables spanning the phase space of the system, in our 
case the Fourier coefficients A 1 -, A 2 -, B 1 -, B 2 , rrij, and Uj, with j — 1, . . . , JT. The condition 
t > ttr allows for the transients to die out. Typically, t tr is rather safely fixed to 1500, which 
correspond to about five years. 

When the Hadley equilibrium is the universal attractor, the volume Vbb is zero, while it 
is non-zero if the computed orbit is attracted to a periodic orbit, a two-torus or a strange 
attractor. In all cases Vbb, which represents the bulk size of the attractor in phase space, 
grows with Te- More precisely, each of the factors in the product (J78)) increases with 
Te, so that expansion occurs in all directions of the phase space. This matches the basic 
expectations on the behavior of a dissipative system having a stronger energy input. 

In the right panel of Fig. ^^we present a plot of log(VBB) as function of Te for the selected 
values of JT = 8, 16, 32, and 64. In the case JT = 8, Vbb obeys with great precision the 
power law Vbb (Te — T E r **) 7 in the whole domain Te > 9. The best estimate for the 
exponent is 7 ~ 40. Given that the total number of Fourier components is 6 x JT = 48, this 
implies that the growth of the each side of the bounding box is on the average proportional 
to about the 5 /6 th power of (Te — T E nt )- 

For higher values of JT, two sharply distinct and well defined power-law regimes occur. 
For JT = 16, in the lower range of (Te — T E nt ) - corresponding in all cases to Te ^ TJP*+1.5 
- the volume of the bounding box increases with about the 35 th power of (Te — T^"*), while 



in the upper range of (Te — T E nt ) - for Te > T E nt + 1.5 - the power-law exponent abruptly 
jumps up to about 80. For JT = 32 the same regimes can be recognized, but the values 
of the best estimates of the exponents are twice as large as what obtained with JT = 16. 
Similarly, for JT = 64 the best estimates of the exponents are twice as large as for JT = 32. 
The results on the power law fits of Vbb oc (Te — T E nt )~< are summarized in Table I1V1 
We emphasize that in all cases the uncertainties on 7, which have been evaluated with a 
standard bootstrap technique, are rather low and total to less than 3% of the best estimate 
of 7. Moreover, the uncertainty of the power-law fit greatly worsens if we detune the value 
of T E rit by as little as 0.3, thus reinforcing the idea that fitting a power law against the 
logarithm of (Te — T E nt ) is a robust choice. 

When considering separately the various sides of the bounding box hyperparallelepiped 
(not shown), i.e., each of the factors in the product (|78jh we have that for JT = 8 all of 
them increase as about (T E - T c E rit f^ in the whole range. For JT = 16, 32, and 64, in 
the lower range of Te each side of the bounding box increases as about the l/3 rd power of 
(Te — Te^), while in the upper range of Te each side of the bounding box increases as about 
the 5 /6 th power of (Te — T E nt ). Selected cases are depicted in the right panel of Fig. EI 
So for a given value of truncation order JT, the ratios between the ranges of the various 
degrees of freedom are essentially unchanged when varying Te, so that the system obeys a 
sort of self-similar scaling with Te- 

Summarizing, for sufficiently high truncation order (JT > 16) a robust parametric de- 
pendence is detected for the volume of the bounding box as a function of Te: 

I 1/3, T E - T c E rit < 1.5, 
V BB oc(TE-T^y, 1 = eN e~ E (79) 

[5/6, T E -T E rU > 1.5, 

where N = Q x JT is the number degrees of freedom. 

The comparison, for, say, JT =16 and 32, of factors in ()78|) having the same order for 
the same value of Te — T E lt provides insight about the sensitivity to model resolution. In 
the following discussion, we examine the variables A 1 - but similar observations apply to all 
other variables A*, Bj, B?, Uj, and rrij. The factors related to the the gravest modes, such 
as [max (A^ (t)) — min (A\ (£))], agree with high precision, thus suggesting that the large 
scale behavior of the system is only slightly affected by variation of model resolution. When 
considering the terms related to the fastest latitudinally varying modes allowed by both 
truncation orders, such as [max (A) (t)) - min (A] (£))] with 17 < j < 32, we have that 
those obtained for JT = 32 are larger than the corresponding factors obtained for JT = 64, 



and the distance between pairs of the same order increases with j. See the right panel of 
Fig. ^2 This is likely to be the effect of spectral aliasing (mJ: the fastest modes of the 
model with lower resolution absorb the dynamics contained in the scales which are instead 
resolved in the higher-resolution model. The same effect is observed when comparing, for 
JT = 16 and 32, coefficients of the same order such as [max (Aj (t)) — min (A* (t))] with 
9 < j < 16. The JT = 8 case does not precisely match this picture. 



IV. STATISTICAL PROPERTIES OF THE TOTAL ENERGY AND LATITUDI- 
NALLY AVERAGED ZONAL WIND 

In this section the model (|61| ) - (j66|) is studied by means of observables (functions of state 
space variables) of physical significance, as opposed to the quantities derived from the Lya- 
punov exponents and the volume of the bounding box used in Sec. IIII CI which are more 
typical indicators used in dynamical systems analysis. 

A. Total energy 

The total energy of the system E(t) defined in ()75|) isa very relevant observable of phys- 
ical significance for the system. In Table |l] we report its conversion factor between the 
non-dimensional and dimensional units. For the Hadley equilibrium, the time-independent 
expression for the total energy is derived by plugging (foT?|) into and then computing the 
integral (|T5|) : 



The total energy is proportional to T| and is mostly stored as potential energy which is 
described by the second term of the sum in (|%U|) . 

In Fig. El we present the results obtained for the various values of JT used in this work. 
In the left panel we present the JT = 64 case, which is representative of what obtained 
also in the other cases. The time-averaged total energy is monotonically increasing with Te, 
but when the system enters the chaotic regime, E(t) is much lower than the value at the 
coexisting Hadley equilibrium. This behavior may be related to the much larger dissipation 
fuelled by the chaos-driven activation of the smaller scales. In the chaotic regime E(t) is 
characterized by temporal variability, which becomes more and more pronounced for larger 
values of T E . 

In the right panel of Fig.^Jwe compare the cases JT = 8, 16, 32 with respect to JT = 64. 
The overall agreement of E(t) is good but progressively worsens when decreasing JT: for 
JT = 32, the maximal fractional difference is less than 0.01, while for JT = 8 it is about 
one order of magnitude larger. Differences between the representations given by the various 
truncations levels also emerge in power law fits such as E(t) oc Tg. In the regime where the 
Hadley equilibrium is attracting, this fit is exact, with exponent 7 = 2. For Te — T E nt < 1.5 
and Te > T^ (the value of the first Hopf bifurcation, see Table ITTT|) . for all the values of 




2 
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JT the power law fit is good, with 7 = 1.9 ± 0.03, so that a weakly subquadratic growth 
is realized. For Te — T E rit > 1.5, only the JT = 32 and 64 simulations of E(t) obey with 
excellent approximation a weaker power law, with 7 = 1.52 ± 0.02 in both cases, while the 
cases JT = 8 and 16 do not satisfactorily fit any power law. 

The agreement worsens in the upper range of Te, which points at the criticality of the 
truncation level when strong forcings are imposed. Nevertheless, the observed differences 
are strikingly small between the cases, say, JT = 8 and JT = 64, with respect to what 
could be guessed by looking at the Lyapunov dimension, entropy production, and bounding 
box volume diagnostics analyzed in the previous sections, where essentially only JT = 32 
and JT = 64 had a satisfactory agreement. This suggests that when analyzing global 
observables, the resolution requirements for obtaining good statistical indicators are much 
more relaxed. 



B. Zonal wind 

We here examine the latitudinal average, denoted by (•), of U and m: 

1 9 ^ TP 

(U(y,t)) = - U(y,t)dy = - — > ( 81 ) 

L In IX * ' 7 

J0 i=i,jodd J 

(m(y,t)) = - m(y,t)dz=- V — . (82) 

L Jo n ■ aa J 

]=1,3 odd 

Since U(y, t) represents the zonal average of the mean of the zonal wind at the two pressure 
levels pi and 793 at latitude y, (U(y,t)) is proportional to the total zonal momentum of the 
atmosphere. Instead, (m(y,t)) represents the spatially averaged halved difference between 
the the zonal wind at the two pressure levels pi and p$. Computation of such space averages 
at the time-independent Hadley equilibrium (|T0|) is straightforward: 

(m (y)) = (U (y)) = ^2 • (83) 

v N 

Since we cannot have net, long-term zonal forces acting on the atmosphere at the surface in- 
terface, the spatial average of the zonal wind at the pressure level p$ must be zero. Therefore, 
the outputs of the numerical integrations must satisfy the following constraint: 



(m (y,t)) = (U(y,t)), 



(84) 



where X denotes the time-average of the field X. The constraint (fMj) is automatically 
satisfied at the Hadley equilibrium. 

The results are presented in Fig. In the left panel we plot the outputs for JT = 64, 
which, similarly to the total energy case, is well representative of all the JT cases. We 
first note that the constraint (J8%j) is obeyed within numerical precision. The average winds 
are monotonically increasing with Te, but, when the system enters the chaotic regimes, the 
averages (m(y, t)) = (U (y, t)) have a much smaller value than at the Hadley equilibrium, and 
they display sublinear growth with T%. Moreover, for Te > T E rit the temporal variability of 
the time series (m(y,t)) and (U(y,t)) increases with Te- The variability of (m(y,t)) results 
to be slightly larger than that of (U(y,t)), probably because the latter is related to a bulk 
mechanical property of the system such as the total zonal momentum. 

Since we are dealing with a quasi-geostrophic system, these observations on the wind fields 
imply that while the time-averaged meridional temperature difference between the northern 
and southern boundary of the system increases monotonically with Te, as to be expected, 
the realized value is greatly reduced by the onset of the chaotic regime with respect to the 
corresponding Hadley equilibrium. This is the s igna ture of the negative feedback due to a 
mechanism similar to the baroclinic adjustment [6i|: when the poleward eddy transport of 
heat is realized, it causes the reduction of the meridional temperature gradient, thus limiting 
the wind shear. Note that in this model the adjustment, as opposed to the general case, 
is essentially correct in a variational context, since only one zonal wave is considered, and 
so the fastest growing, unstable wave is also the wave transporting northward the largest 
amount of heat 0,0] . Nevertheless, the adjustment mechanism does not keep the system 
close to marginal stability, as envisioned in some baroclinic adjustment theories, since for 
Te > T E rit both the instantaneous and the time-averaged fields of the system are completely 
different from those realized at the Hadley equilibrium. 

The effects of lowering JT are illustrated in Fig. EH right. The overall agreement, ex- 
pressed by a small value of the fractional differences, progressively worsens for smaller JT. 
Notice the similarity of the functional shapes with Fig. II 21 right. The results in Fig. IT3*1 right 
can be summarized as follows: the coarser-resolution models have higher total temperature 
difference between the two boundaries for values of Te up to about 30 and lower temperature 
differences for higher values of Te- This implies that while for Te < 30 the latitudinal heat 
transport increases with JT as a positive trade-off between the higher number of unstable 
baroclinic modes (within a sloppy linear thinking) or, better, smaller scale baroclinic con- 
version processes taking place in a higher-dimensional attractor, and the enhancement of 
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Again, differences between the various truncations levels emerge as one attempts power 
law fits of the form (m(y,t)) = (U(y,t)) oc T^. For the Hadley equilibrium regime we have 
7 = 1. For Te < 10 and above the first Hopf bifurcation, for all values of JT the power law 
fit is good, with 7 = 0.875 ± 0.05. For Te — T E nt > 1.5, only the simulations with JT = 
32 and 64 obey a power law (with 7 = 0.58 ± 0.02) with excellent approximation, while the 
realizations of the JT = 8 and 16 cases do not fit any power law. 

By examining more detailed diagnostics on the winds, such as the time-averaged lati- 
tudinal profiles of U(y) and of m(y) (Fig. 1X4^ . relevant differences are observed between 
JT = 8 and the other three cases. Results are presented for JT = 8 and JT = 32, the 
latter being representative also of JT = 16 and 64. We first note that already for Te = 
9 and 10, such that only a weakly chaotic motion is realized, the U(y) and m(y) profiles 
feature in both resolutions relevant qualitative differences with respect to the corresponding 
Hadley equilibrium profile, although symmetry with respect to the center of the channel is 
obeyed. The U(y) and m(y) profiles are different (the constraint (J84|) being still satisfied), 
with U(y) > m(y) at the center and U(y) < m(y) at the boundaries of the channel. Never- 
theless, like for the Hadley equilibrium, both U(y) and m(y) are positive and are larger at 
the center of the channel than at the boundaries. Consequently, at pressure level p\ there is 
a westerly flow at the center of the channel and easterly flows at the two boundaries, and 
that at pressure level P3 the wind is everywhere westerly and peaks at the center of the 
channel. Such features are more pronounced for the JT = 32 case, where the mechanism of 
the convergence of zonal momentum is more accurately represented. 

For larger values of Te, the differences between the two truncation levels become more 
apparent. For JT = 8, the observed U(y) and m(y) profiles tend to flatten in the center 
of the channel and to become more similar to each other. Therefore, somewhat similarly 
to the Hadley equilibrium case, the winds at the pressure level pi tend to vanish and all 
the dynamics is restricted to the pressure level p 3 . The m(y) profiles for JT = 32 are quite 
similar to those of JT = 8, even if they peak and reach higher values in the center of the 
channel and are somewhat smaller at the boundaries. So when a finer resolution is used, a 
stronger temperature gradient is realized in the channel center. The U(y) profiles obtained 
for JT = 32 are instead very different. They feature a strong, well-defined peak in the 
channel center and negative values near the boundaries. Therefore, the winds in the upper 
pressure level are strong westerlies, and peak in the center of the channel, while the winds in 
the lower pressure level feature a relatively strong westerly jet in the center of the channel 
and two compensating easterly jets at the boundaries. The fact that for higher resolution 
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more efficient mechanism of barotropic stabilization, which, through zonal wind convergence, 



Examination of the latitudinal profiles in Fig. El clarifies our choice to extend the latitudi- 
nal domain of the model beyond the geometrically and geographically realistic mid-latitude 
channel. Thanks to this, the wind fields in the central portion of the domain (the latter 
corresponds to mid-latitudes and is of primary interest in this work), are rather different 
than at the boundary regions. The observed features, and especially the presence of a jet, 
are in qualitative agreement with the real atmosphere if models having truncation order of 
JT > 16 are used. 

Summarizing, by considering the latitudinal average of the wind fields in the mid-latitudes 
range [0.251^,0.751^], the JT = 8 model greatly differs from the higher resolution models, 
since (m) and especially (U) are underestimated. Indeed, these diagnostics do not only rely 
on a global balance, which is relatively weakly resolution-dependent (see previous section), 
but also on the resolution-sensitive representation of internal processes such as the zonal 
wind convergence. 



keeps the jet together [37|, |38j, |39 




V. SUMMARY AND CONCLUSIONS 



We have described the construction and the dynamical behavior of an intermediate com- 
plexity model of the atmospheric system. The ab-initio equations of dynamics and thermo- 
dynamics of a stratified fluid are specialized to the quasi-geostrophic motion and a new de- 
tailed derivation of the quasi-geostrophic two-layer model of the planetary scale atmospheric 
flow in a mid-latitudes beta-plane is provided. The derivation is performed by retaining, 
at each step, the variables as expressed in physical units, while the non-dimensionalization 
procedure, useful for the numerical integrations, is introduced at last. 

A single zonal wave solution is assumed and a partial differential equation is derived 
for its coefficients. By a spectral discretization in the latitudinal direction (using a Fourier 
half-sine expansion), the latter equation is reduced to a system of N = 6 x JT ordinary 
differential equations, where JT + 1 is the number of nodes of the (latitudinally speaking) 
fastest varying base function. We have considered the cases JT = 8, 16, 32, and 64. 

By increasing the parameter Te, corresponding to the imposed equator-to-pole tempera- 
ture gradient, the system develops a strange attractor in phase space. The route leading to 
the formation of this strange attractor involves: 

• a Hopf bifurcation at Te = T§ responsible for the loss of stability of the Hadley 
equilibrium (corresponding to corresponding to baroclinic instability), where a periodic 
orbit branches off; 

• a Hopf bifurcation where a two-torus is created; 

• a finite number of quasi-periodic period doublings of the invariant two-torus; 

• two-torus breakdown at Te = T E rit . 

Statistical indicators, such as lagged autocorrelations, have been used to characterize the 
observed quasi-periodic or strange attractors for various values of Te- To generate the 
required time series, a physically relevant observable has been computed, the total energy of 
the system. For Te close to T E rit quasi-periodic intermittency and very weak chaoticity are 
detected. The corresponding flow pattern might be classified as an amplitude vacillation, 
like for the two-torus dynamics. For larger Te the lagged autocorrelation typically decays 
(exponentially) fast. The observed route to chaos is qualitatively the same for JT = 16, 32, 
and 64, and the values T E and TJT 1 * weakly depend on JT (Table ITTTj) . Structural differences 
occur for JT = 8: a transition to a three-torus, yielding a modulated amplitude vacillation, 



The strange attractor is further studied by means of the Lyapunov exponents, where we 
have varied both Te and model resolution JT. Although the system qualitative behavior 
is analogous for different values of JT, there are significant quantitative differences. In all 
cases, the maximal Lyapunov exponent Ai increases with Te, and it is possible to robustly fit 
a power law of the form Ai oc (T# — Tjr 1 *) 7 . For Te fixed, the maximal Lyapunov exponent 
decreases with JT (so that the predictability time increases). On the contrary, the metric 
entropy increases linearly with Te for all examined values of JT, and is larger for larger 
values of JT. In other words, the fastest (the total) dynamical instability of the system is 
smaller (larger) for larger JT, where the dynamics is more accurately represented. 

The Lyapunov dimension increases with both Te and JT. The dependence of on 
Te is qualitatively the same for all values of JT: by increasing Te there is an initial phase 
where the dimension quickly grows with a power law Dl oc (Te — T E nt y< , followed by a 
linear scaling regime. For large Te, the dimension saturates and depends sublinearly on Te- 
The latter effect is, of course, more evident for small values of JT. It provides a measure of 
accuracy of the spectral discretization (as far as the details of the dynamics are concerned), 
which turns out to depend on Te- 

When considering the bounding box of the system, i. e. the minimal hyperparallelepiped 
containing the attractor in phase space, for sufficiently high truncation order JT each side 
of the box increases as oc (T E - Tjf*) 1 / 3 for T E - Tjf* < 1.5 and as oc (T E - T^"*) 5 / 6 for 
larger values of Te- So for a given value of JT the ratios of the ranges of the various degrees 
of freedom remain essentially unchanged when varying Te, yielding a self-similar scaling 
property. The volume of the bounding box Vbb then results to increase as oc (Te — T E nt ) N ^ 3 
and as oc (T# — j^™*) 5 ^ 6 j n the mentioned domains of Te- 

A peculiar feature of this dynamical system is the rather smooth dependence on the pa- 
rameter Te of all the examined properties of the strange attractor. No windows of periodicity 
have been detected in the chaotic range and this is quite uncommon especially when com- 
paring with low-dimensional chaoti c sy stems such as the Henon-Pomeau mapping jsil 0] 



or the Lorenz flow 



[loj l (also see 0, 



491]). Althou gh s tructural stability 
question, other stability concepts (such as robustness |55j) might provide an alternative and 
more practical theoretical basis for the explanation of the observed parametric smoothness, 
perhaps also for other systems of intermediate and high dimensionality. 

Despite the sensitivity of Lyapunov exponents and dimension to model resolution JT, 
certain observables of physical interest, such as the time-averaged total energy of the system, 
or the time-averaged spatially averaged zonal wind fields, are in quantitative agreement for 
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global balances, which turn out to be only slightly affected by model resolution. When the 
system enters the chaotic regime, the average total energy and average zonal winds have 
lower values than those of the coexisting - and unstable - Hadley equilibrium, because the 
chaos-driven occupation of the faster-varying latitudinal modes fuels viscous dissipation, 
which acts preferentially on the small scales. Other mechanisms which are present in the 
real atmosphere, such as the barotropic governor 65| . are not represented in this schematic 
model. Moreover, the total energy and the average wind field at the Hadley equilibrium 
depend quadratically and linearly on T E , respectively, in the chaotic regimes such quantities 
obey a subquadratic and sublinear power law oc T^, respectively. For both quantities, the 
exponents of the power laws decrease abruptly as T# — T| ot crosses 1.5. An analogous sharp 
change is observed for Vbb, which suggests the onset of a self-similar scaling law. 

Nevertheless, when analyzing more detailed diagnostics on the winds at the two pressure 
levels, relevant differences emerge between the model with JT = 8 and those with the 
higher resolutions. For JT = 8 the wind profiles are rather flat in latitude and very weak in 
the lower pressure level. For the higher resolution models the winds in the upper pressure 
level are strong westerlies and peak in the center of the channel (corresponding to mid- 
latitudes), in qualitative agreement with reality. The winds in the lower pressure level 
feature a relatively strong westerly jet in the center of the channel and two compensating 
easterly jets at the boundaries. The fact that for higher resolution the wind profiles are less 
smooth and have more clear-cut jet-like features is related to the more efficient mechanism 
of barotropic stabilization, which, through zonal wind convergence, keeps the jet together. 

The model we study, although admittedly very schematic, is Earth-like in that it features 
some fundamental processes determining the general circulation of the Earth atmosphere, 
in particular: 

• the complex process of atmospheric baroclinic conversion, transforming available po- 
tential energy associated with (latitudinally) differential Sun heating into kinetic en- 
ergy of synoptic scale motions of the mid-latitudes atmosphere; 

• nonlinear stabilization by eddy momentum convergence from non-symmetric baroclinic 
disturbances into the zonal jet; 

• viscous dissipation. 

While the baroclinic conversion process is essentially well represented in all models (even if 
those with higher resolution are more efficient in the conversion for large Tg, since conver- 
sion can take place also on smaller spatial scales), the descriptions of the barotropic zonal 



wind convergence and of the viscous dissipation are much more critically dependent on the 
latitudinal truncation order, since the latter processes are represented by terms involving 
the latitudinal derivatives of the fields. 

When a larger pool of available energy is provided, the dynamics of the system is richer, 
since the baroclinic conversion process can transfer larger amounts of energy to the distur- 
bances: for each given value of JT, the largest Lyapunov exponent, the metric entropy, the 
Lyapunov dimension of the attractor, the mean and the variability of the total energy and 
of the latitudinally averaged zonal wind fields all increase with T#. The enhancement of the 
efficacy of the baroclinic conversion process when higher resolution is adopted is highlighted 
by the increase with JT, for a fixed value of Tg, of the number of linearly unstable modes 
of the Hadley equilibrium, of the Lyapunov dimension of the attractor, and of the metric 
entropy. 

The critical dependence of the efficiency of the two mentioned stabilizing processes on the 
model resolution is illustrated by several results, e.g. in the dependence of the parameters 
Tfi and T E rtt on JT (it is easier to destabilize a system with lower resolution), in the fact 
that there are fewer unstable modes of the Hadley equilibrium for larger values of JT in the 
vicinity of Tjf, in the fact that the predictability time increases with JT for a given value 
of Te and in the features of the latitudinal profiles of the winds. 

Although relevant ingredients of geometrical (horizontal convergence due to the Earth 
curvature, latitudinal boundary conditions at the margins of the middle latitude circumpolar 
vortex, etc.) and dynamical (stabilization mechanisms other then momentum convergence 
such as the so-called barotropic governor 0|) nature of the real atmospheric circulation are 
still missing in this, very preliminary, theoretical representation, some important general 
conclusions are drawn from the described results. 

Pessimistic conclusions (in increasing order of pessimism): 

• No simple mean field or macroscopic adjustment theory can be formulated for such 
complex nonlinear systems, even for relatively simple models as those proposed in this 
work. 

• It is, in general, doubtful whether invariant manifolds in phase space - such as fixed 
points, periodic orbits - carry any useful information concerning the general circulation 
of the system. 

• Beyond the time of deterministic predictability, "averaging" is of no practical use; it 
is not clear what else should be done in order to produce useful - in a statistical sense 
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Optimistic conclusions: 

• Although some dynamical system properties, such as Lyapunov exponents and di- 
mension, are strongly model-dependent, some other - of great physical interest - are 
not. 

• Increasing refinement (number of degrees of freedom) of models may produce smoother 
dependence on macroscopic parameters. 

• It is not outside the range of practically feasible, although possibly challenging, 
projects to put together an intermediate dimensionality model - with hundreds of (well 
chosen!) degrees of freedom - with stable properties which is relevant for a theory of 
general atmospheric circulation. 
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APPENDIX A: ON THE NUMERICAL METHODS 



We begin by describing the projection operator used in the definition of the vector 
field (pjTjl - tpjj) . As it is customary with climatological s pec tral models 0], a pseudospectral 
method is used, also known as Fourier collocation 

The fields A, B, U, m, appearing in the nonlinear terms of (}6*T]) -()66 j) . are first evaluated 
at JT collocation points yi, . . . , jjjt, equally spaced in the ^/-domain (0, L y ). This is achieved 
by a Discrete Sine Transform of Aj, A 2 , Bj, B 2 , Uj, rrij, with j = 1, . . . , JT. The terms 
involving second derivatives with respect to y are also computed in this way, by premulti- 
plying for a suitable coefficient involving the wave numbers Wj. Then all the nonlinear terms 
are evaluated pointwise, at each of the collocation points yi, . . . , yjT- Lastly, an inverse Dis- 
crete Sine Transform is carried out, yielding the Fourier coefficients of the nonlinear terms. 
The software library fftw3 publicly available at www.fftw.org, has been used for the 
Discrete Sine Transform. 

The numerical solution of the system of ordinary differential equations (}6Tj) - (j66j) is com- 
puted by means of a standard Runge-Kutta-Fehlberg(4,5) algorithm 6|J with adaptive step- 
size, where the approximated solution is carried by the order five method. The local trun- 
cation error is kept below l.e — 6. The stepsize adjustment procedure is similar to that of 
D0PRI5, available at (www.unige.ch/~hairer). 

The total energy of (}6*T]) -()66 p . is computed according to (|75|). In terms of the Fourier 
coefficients Aj, B{, . . . , this yields the expression 
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For the computation of the averages in Sec. IIVI time series of 315360 adimensional time 
units (1000 years in natural units) have been computed for all values of JT, preceded by 
a transient of five years (time is expressed in the scale of the system, see Table UJ). The 
observables E(t), U(y,t), and m(y,t) have been sampled every 0.216 time units (four times 
a day), thereby obtaining time series of 1460000 elements. The sample mean and sample 



standard deviation have been computed according to the usual formulas: 

E^) = -^E(U), 4 = 7 $> 2 (t 4 )-^ 2 . (A2) 
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The initial condition used for all computations is A\ = —0.8, A\ = 0.65, B\ = 0.2, B\ = 
0.2 , = 0.4, .Bf = 0.1, Z7i = 1.26, mi = 1.1, with the remaining coefficients set to 0, as 



in 



24|. 



APPENDIX B: LYAPUNOV EXPONENTS 



The Lyapunov exponents of system (|6i p -([66 |) are estimated according to the algorithm 
described by Galgani, Giorgilli, Benettin and Strelcyn 7(|. The first variational equations 
of are integrated during a period of time T, with the identity matrix as initial 

condition. During integration, at time t the canonical orthonormal basis is mapped onto a 
new set of vectors (v\, w\, ■ ■ ■ , v^r), where N = 6 x JT is the dimension of the phase space. 
Each vector tends to align itself along the direction of maximal expansion (or of minimal 
compression). Thus all v*'s tend to collapse onto one direction. To prevent this, the Gram- 
Schmidt process is applied to (v* 1 , v^ 1 , . . . , v^) at t = t\, yielding a set (vj 1 , . . . , v^) of 
orthogonal vectors. The vectors are normalized by putting w* 1 = v*- 1 / 1 1 v- 1 1| for j = 1, . . . , N. 
Then a new frame of vectors (v^ 2 , . . . , v^), with t<i = 2t±, is computed by integrating the first 
variational equations taking as initial condition the orthonormal vectors (w|, . . . , wjy) from 
the previous step, and the whole process is repeated. At iteration step k, define = kt\ 
and 

k 

\ ||v* fe || and w* fc = — f— for j = 1, . . . , N. 



c k 

7=i ' ' H v i 



The orthonormalization process does not change the direction of v| fe , so that still points 
to the direction of maximal stretch. Denoting by Xj, j — 1, . . . , N, the Lyapunov exponents 
in decreasing order of magnitude, the length c\ of v{ k is approximately proportional to e kXl . 
The plane spanned by v{ k and v 2 fe is not changed by the Gram-Schmidt process and tends 
to adjust to the subspace of maximal growth of surfaces. The rate of growth of areas is 
proportional to e k<yXl+X2 \ In particular, since v' fc = w^ fe and are orthonormal, the length 
of the projection of v 2 fe upon w^' is proportional to e kXi . A similar argument for growth 



n k\i 

o piupui Liunai \jU t 

estimated by the averages 

A, « 

k 

where k = T jt\. We have chosen T of the order of 3150 adimensional time units (about 10 



of volumes yields that c k is proportional to e fcAj . Therefore, the Lyapunov exponent is 

Xj* -log^), (Bl) 



years in natural units) for all values of JT, while t\ has been chosen as 0.864 adimensional 
time units (1 day), which allow for an excellent convergence of the exponents. 

Actually, we have used a version of the algorithm 70j in which the variational equations 
are not integrated explicitly, but approximated by means of numerical differentiation: N 
trajectories are simultaneously integrated, starting from points nearby a reference orbit. 
The distances from the reference orbit are normalized at regular time steps \(\ - 

The library LAPACK fwww.netlib.orsr") has been used for Gram-Schmidt orthoeonaliza- 



tion and for other computations in this work. 
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TABLE I: Variables of the system and non-dimensionalization factors. For A\, B 2 , m< 

U n , and w n , the index n ranges from 1 to JT. For Xj, n = 1, . . . , 6 X JT. 
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TABLE II: Values of the parameters used in this work and non-dimensionalization factors. 



JT 


T# 


rpcrit 
1 E 
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7.83 


9.148 


16 


8.08 


8.415 


32 


8.28 


8.522 


64 


8.51 


8.663 



TABLE III: Approximate values of the parameter Te where the Hadley equilibrium loses stability 
via Hopf bifurcation (T E ) and where the onset of the chaotic regime occurs (T E rit ) for each of the 
considered orders of truncation JT. See text for details. 



JT 


j[loi 




< 0.5] 


7 [log(T B - T E rU ) > 0.5] 
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40 ± 1 




40 ± 1 


16 




33 ±3 
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320 ± 1 



TABLE IV: Powerdaw fits of the volume of the bounding box as Vbb oc (Te — T E " lt )' y in two 
different ranges of Te — T E lt for each of the considered orders of truncation JT. See text and 
Fig- EH for details. 
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FIG. 1: Sketch of the actual geographical area corresponding to the simplified (3 channel. The 
local x and y directions and the /3-channel width L y are indicated. The mid-latitudes range from 
l/4Ly to 3/ALy, corresponding to a 45° latitudinal belt centered at 45° N. 
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FIG. 2: Sketch of the vertical-longitudinal section of the system domain. The domain is periodic 
in the zonal direction x with wavelength L x . At each pressure level, the relevant variables are 
indicated. 
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FIG. 3: Number of linearly unstable modes at the Hadley equilibrium as a function of the param- 
eter T E for JT = 8, 16, 32, 64. 
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FIG. 4: Time-evolution of the component A\ of system (j61[) - ()66[) . starting from the initial condition 
mentioned in Appendix Left: Te = 8.484631. Right: Te = 8.484632. No transient has been 
discarded. 



FIG. 5: Left: projection on of an orbit on the attracting two-torus of l|61 [) -(|66 [) for 

Te = 8.484632. Right: same as Left, projection on (A\,Ul). Units as indicated in Table [I] A 
five-year transient has been discarded. The phase-space region where the orbit accumulates more 
densely is due to intermittency of saddle-node type near the location of the periodic orbit occurring 
for Te = 8.484631, compare Fig. 0]and see text for details. 




FIG. 6: Projections on of a Poincare section of the attractor of (|61 [) -(|66 [) . obtained 

by intersecting it with a hyperplane Ui = cq for several values of Te- From (A) to (F) Te 
is, respectively, 8.516, 8.52, 8.521, 8.522, 8.58, 10. The value Co of the section is fixed at 0.66 
(A) to (D) and is 0.7 and 0.8 for (E) and (F) respectively. Also notice the different axis scale for 
the last two plots. 
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FIG. 7: Autocorrelations of the total energy time series on the attractor of (|61 )) -(|66 )) . for various 
values of T E . Left: T E = 8.521, 8.522, 8.58; Right: T E = 9, 10, 18. 




FIG. 8: Left: Lyapunov exponents Xj for JT = 32 and for j = 1, 2, 32, 96, 160, 191, 192 as a 
function of T E . Right: Spectrum of the Lyapunov exponents for T E = 9, T E = 30, and T E = 110. 
Units as for Xj and T E as described in Table [I] 
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FIG. 9: Lyapunov dimension of the attractor of (|61|) - (|66j) as a function of Te for JT = 8, 16, 32, 
and 64. All the straight lines are parallel and the domain of validity of the linear fit is apparently 
homothetic. 
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FIG. 10: Left: maximal Lyapunov exponent on the attractor of (|6H) - (j66j) as a function of Te for 
JT = 8, 16, 32, 64. Center: Log-Log plot of the predictability time of the system t p = X^ 1 versus 
Te — T E nt . Power laws (t p oc (Te — T^p*) 7 ) are detected for all considered values of JT. Right: 
metric entropy. Linear dependences h ~ (3(Te — Tg lt ) occur for all values of JT. 




FIG. 11: Left: Volume of the bounding box Vg^ of the attr actor as a function of the detuning 
parameter Te — TJP* for JT = 8, 16, 32, 64. For description of the power law fits, see text and 
Table Hvl Right: Value of the corresponding sides of the bounding box pertaining to the variables 
Aj for JT = 32 (red lines) and 64 (magenta lines). Notice the two power-law regimes mentioned 
in the text. 




FIG. 12: Left: E(t) for the Hadley equilibrium (black line) and deduced from the observed fields 
in the chaotic regime for JT = 64 (magenta line) ; the magenta dashed line delimit the ex-confidence 
interval. Right: fractional deviations of E(t), for JT = 8, 16, and 32, with respect to JT = 64. 
See text for details. 




FIG. 13: Left: (U) = (m) for the Hadley equilibrium (black line) and deduced from the observed 
fields in the chaotic regime for JT = 64 (magenta line); the magenta dashed and dotted lines 
delimit the cr-confidence interval for (U) and (m), respectively. Right: fractional deviations of 
(U) = (m) for JT = 8, 16, and 32 with respect to JT = 64. See text for details. 
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FIG. 14: Time-averaged latitudinal profiles U(y) (solid lines) and m(y) (dashed lines). In all 
figures the black solid line indicates the U(y) = m{y) profile of the Hadley equilibrium, the blue 
and red lines refer to the cases JT = 8 and JT = 32, respectively. The values of are indicated. 
Note that the vertical scale for Te = 9 and 10 is about 1/2 as for the other two figures. 
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